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Summary 

This report reviews two topics that are important to discrete element method (DEM) modeling the 
behavior of soils (such as lunar soils): (1) methods of modeling particle shapes and (2) analytical 
representations of particle size distribution. These two topics are somewhat unrelated, but both are 
reviewed in this report as a part of ongoing efforts to develop models of lunar soil mechanics. The choice 
of particle shape complexity is driven primarily by opposing tradeoffs with total number of particles, 
computer memory, and total simulation computer processing time. Also, the choice is dependent on 
available DEM software capabilities. For example, PFC2D/PFC3D and EDEM support clustering of 
spheres; MIMES incorporates superquadric particle shapes; and BLOKS3D provides polyhedra shapes. 
Most commercial and custom DEM software supports some type of complex particle shape beyond the 
standard sphere. Convex polyhedra, clusters of spheres and single parametric particle shapes such as the 
ellipsoid, poly-ellipsoid, and superquadric, are all motivated by the desire to introduce asymmetry into the 
particle shape, as well as edges and corners, in order to better simulate actual granular particle shapes and 
behavior. 

An empirical particle size distribution (PSD) formula is shown to fit the example desert sand data 
from Bagnold, using a coarse-grade exponent c and small-grade exponent 5. An additional empirical 
parameter y has been introduced to connect the two (log-log) linear fits in the small- and coarse-grade 
particle size regions. Note that the PSD has been defined so that the integral from size 0 to °° is equal to 
100 in units of percent. Particle size data of JSC- la obtained from a fine particle analyzer (FPA) at NASA 
Kennedy Space Center is also fitted to a similar empirical PSD function in order to extract the equivalent 
small- and coarse-grade exponents. 


1.0 Introduction 

The starting point with any discrete element method (DEM) software is to model grains as spheres of 
a specified diameter range. The spherical particle shape greatly simplifies simulations and accommodates 
the maximum number of particles for any given central processing unit (CPU) execution time budget. 
However, it is well known that simulations of most real granular assemblies, such as sand or lunar 
regolith, are not accurately modeled by spheres because of the rolling of contact surfaces versus the actual 
frictional sliding of real particle surfaces. Most researchers active in granular modeling agree that a 
particle shape more complex than the sphere is required to generate simulations that can match the 
observable properties of real granular systems (Refs. 1 to 4). Section 2.0 of this report addresses particle 
shape. 
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The particle size distribution is one of the primary factors affecting how a soil behaves mechanically. 
In order to model soil mechanics, it is important to start out with an understanding of the size distribution 
of the soil that is being modeled. Section 3.0 of this report examines analytical models that can 
approximate a lunar soil particle size distribution. Appendix A shows magnified images and elemental 
spectra of JSC- la lunar regolith stimulant, developed by the NASA Johnson Space Center. It should be 
noted that particle size distributions have been measured by several methods, which do not always agree 
with one another because each method is measuring a different aspect of particle “size.” Therefore, these 
analytical functions should not necessarily be considered definitive for lunar soil or lunar simulants until 
further work is accomplished. Instead, these should be considered as tentative “models” of the lunar soil 
particle size distribution that will be refined with further research. 


2.0 Methods of Modeling Particle Shapes 

Here, Subsections 2.1 to 2.7 cover different particle constructs. Subsection 2.8 describes a contact 
detection scheme, given nonspherical particles. Appendix B presents a summary of particle shape types 
associated with various DEM software packages, comparing some of the pros and cons of each of these 
idealized particle characteristics. 


2.1 Clusters of Spheres 

Popular DEM software, such as that developed and marketed by the companies Itasca and DEM 
Solutions Ltd., make use of sphere clusters (Fig. 1). Normally, the spheres making up a sphere cluster 
form a rigid body where the adhesion force between spheres is set to infinity. In some cases, researchers 
prefer to connect the spheres using finite adhesion forces in order to model crushing and breaking of large 
complex granular particles. 

Another variation of particles composed of sphere clusters uses overlapping spheres (Fig. 2) (Refs. 7 
to 9). Even though this is the most general case of sphere clusters, it is not necessarily the most common 
approach because of complications caused by its implementation in most of the commercial DEM 
software. DEM software typically models particles using a constant mass density. If particles of equal 
density overlap, the density doubles. Therefore, it is often necessary to implement a careful scaling of 
mass density of each of the component spheres so that the final cluster density uniformity is optimized. 

An algorithm to adjust mass density in two dimensions is described by Ashmawy (Ref. 8). 



Figure 1 . — Particle shape defined by cluster of 
nonoverlapping spheres. 


Figure 2. — Particle shape defined by cluster of 
overlapping spheres. 
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The equation of the surface a single sphere of the cluster contains in general, four parameters, the 
radius R k and the vector offset r k = (x k ,y k ,z k ) relative to a local particle origin: 
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For a cluster of N spheres, there will be 1 + 4 (N- 1) parameters since one of the spheres is defined to be 
centered at the origin of the local particle coordinate system. If the specific DEM software requires a mass 
density scaling of the component spheres, an additional parameter per sphere is needed, for a total of 
5 N- 3 parameters. If however, a density adjustment algorithm is utilized, such as that proposed by 
Ashmawy, only the final density of the cluster is required for a total number of parameters equal to 4 N - 2. 

Equation (1) is valid for both overlapping and nonoverlapping clusters. A more useful representation 
of the surface equation is given in spherical coordinates as 

x(0, cp) = R k sin 0 cos cp + x k 

y(0, cp) = R k sin 0 sin cp + y k (2) 

z(0, cp) — Rfr cos 0 + z for 0 < cp < 271, O<0<7i 

2.2 Ellipsoid 

Other popular particle types are the two-dimensional ellipse and three-dimensional ellipsoid (Fig. 3) 
(Ref. 10). According to many researchers, the asymmetry provided by the ellipsoid yields the best 
tradeoff between particle complexity and total simulation execution time. The equation for the ellipsoid 
surface is similar to Equation (1): 



where a , b , and c are the half- widths along each of the principal axes. The parametric form of 
Equation (3) is 

x(0,cp) = tfsinOcoscp 

y(0, cp) = b sin 0 sin cp (4) 

z(0, cp) = c cos 0 for 0 < cp < 271, 0 < 0 < n 



i o oj co -is 


X 

Figure 3. — Particle shape defined by an ellipsoid 
(Eq. (3)), with a = 1 , b = 0.5, and c = 2. 
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Neither the ellipse or ellipsoid lends itself to clustering strategies, since the attraction to clustering is 
based on a simple unit particle — the sphere. If the unit particle is an ellipse, the difficulties that arise from 
the ellipse tend to negate the advantages of clustering. 

The ellipsoid and sphere are subsets of the superquadric shapes (Ref. 11). 


2.3 Superquadric 

The equation for the superquadric surface is similar to Equation (3): 
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( 5 ) 


where again a , b , and c are the half-widths along each of the principal axes (Fig. 4). The parametric form 
of Equation (3) is: 


x(0, d) = a |sin 0 cos §\ 2/K sgn(sin 0 cos (]>) 
y (0, d) = b |sin 0 sin §\ 2IL sgn(sin 0 sin (j>) 

( 6 ) 

z(0, d) = c |cos 0| 2/M sgn(cos 0) 

for 0<d<27t, O<0<71 

where sgn(x) is the sign of the argument x. The exponent parameters K , Z, and M adjust the blockiness of 
the shape along the three principal axis. If the exponent is less than 2, the shape becomes concave along 
that axis; it becomes convex for an exponent greater than or equal to 2. The superquadrics provides an 
additional degree of asymmetry over the ellipsoid. The greatest degree of asymmetry is achieved by 
defining multiple superquadric parameters in different regions. This strategy leads to the asymmetric 
superquadric. 



Figure 4. — Superquadric particle shape defined 
by Equation (5) where a=1,b=1,c=1, 

K = 0.6, L = 4, and M = 2. 
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2.4 Asymmetric Superquadric 

When Equation (6) is used with P sets of parameters in P regions, an arbitrarily asymmetric shape 
be generated (Fig. 5): 


.2! K\ 


a x |sin 0 cos (|)| sgn(sin 0 cos (|)) 


x(Q,$) = < 


2/K P 

a P |sin 0 cos (])| sgn(sin 0 cos b) 

21 L\ 

b x |sin 0 sin (|)| sgn(sin 0 sin b) 


for 0<([) <()>!, O<0<0! 
for §p_i < b < 271, 0p_! < 0 < 71 
for 0 < (|) < 4>i , O<0<0! 
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2/L P 

b p |sin 0 sin c|)| sgn(sin 0 sin c|>) 

. 2 1 M \ 

c x |cos0| sgn(cos0) for 


for b p_i < (|) < 271, 0p_! < 0 < 71 


o < c[) < 4>i 9 o<0<0! 


z(0,())) = 


2/m p 

Cp |cos0| 


sgn(cos0) 


for § P _i < (|) < 271, 0 P _i < 0 < 7i 


For example, Figure 5 is generated by defining P = 2 regions: 

Region 1: {o < b < 27i, O<0<^7i} 

Region 2: {o < < 2n, ±n < 0 < 71 } 

In the general case, the slopes are not continuous and are characterized by a sharp surface ridge. 



Figure 5. — Asymmetric superquadric particle shape 
defined by Equation (7) where ai = a 2 = 0.7, 
hi = b 2 = 1 , ci = 0.75, c 2 = 2, Ki = K 2 = 4, 

L\- L 2 - 2, Mi = 4, and M 2 = 0.7. 
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2.5 Eight-Quadrant Superquadric 

Defining the eight quadrants of the asymmetric superquadric as three-dimensional regions generates a 
useful particle shape definition: 


*(6,<|>) = 


X9,<l>) = 


z(0,^)) = 


J< 2 1 (sin 0 cos(|)) 2/ ^ 1 

|a 2 |sin0cos(|)| 2/i:2 

[b x (sin0sin(|)) 2/Zl 
[z? 2 |sin0sin <\>\ 2/Ll 

\c x (cos0) 2/Ml 
|c 2 |cos0| 2/M2 


for --|<<|><-| 
for y<(|)<^ 

for 0 < (f) < 7T 
for 7T < 4> < 27T 

for O<0<y 
for y<0<7l 


O<0<71 


O<0<71 


( 8 ) 


The description of the eight-quadrant superquadric form Equation (8) is implemented in Mathematica 
(Ref. 12) and is shown in Figure 6, as well as in Appendix C where the eight-quadrant superquadric shape 
parameters are randomly modified to demonstrate the corresponding range of particle shapes. 
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Figure 6. — Eight-quadrant superquadric video sequence using Equation (8). 
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Figure 7. — Poly-ellipsoid particle shape defined 
by Equation (9) where ai = 2, a 2 = 0.7, hi = 0.5, 
b 2 = 1, ci = 0.75, c 2 = 2, and Ki = K 2 = Li = /_ 2 
= = M 2 = 2. 


2.6 Poly-Ellipsoid 

Poly-ellipsoids, as defined by Peters (Ref. 13), can be viewed as a special case of the eight-quadrant 
superquadric discussed above. In this case (Fig. 7), there are again eight distinct regions corresponding to 
the eight spatial quadrants: 


for - f < (|) < f 
for f < (|) < -y- 

for 0 < c|) < 7T 

(9) 

for 7i < c|) < 2 tt 

for 0<9<f 
for y<0<7l 

where, from Equation (8), L\= L 2 = M x = M 2 = N\ = N 2 are set equal to 2 only in the special case of the 
poly-ellipsoid. This special case defines the poly-ellipsoid and guarantees that the slopes at all quadrant 
boundaries are continuous. 

Care must be taken in defining the parameter sets so that the surface is continuous and does not 
contain voids. A logical method to ensure surface continuity is accomplished by use of the poly-ellipsoid 
particle shape. 



a x sin 0 cos 
a 2 sin 0 cos 

b x sin 0 sin (|) 
b 2 sin 0 sin (|) 

C‘i COS^) 

C 2 COS<\) 


2.7 Polyhedra 

All particle shapes up to this point have been described by closed-form expressions, which is a great 
convenience for computing first order (mass and volume) and second-order (moments of inertia) particle 
parameters of interest for DEM simulations, as well as computing particle interaction (contact forces). All 
simulated particle shapes from compound sphere clusters to poly-ellipses have one goal in common: to 
include shape asymmetry and sharper edges and corners for the purpose of modeling more realistic 
particle shapes. 
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A convex polyhedron can be defined by V vertices, E edges, and F faces, described conveniently by 
the Euler characteristic 


V-E + F = 2 (10) 

The most convenient and compact method to define a general convex polyhedron is by a linear 
system of inequalities, which can be expressed in matrix form as follows: 

A r <b (11) 

where r is an TV x 1 position vector in TV-dimensional space; A is an F x TV matrix (F is the number of 
polyhedron faces); and b is an TV x 1 vector. 

In three-dimensional space r = (x y z ) T , where T denotes “transpose.” The inequality of 
Equation (11) discriminates between points in the polyhedron from those outside of the surface, which are 
not part of the polyhedron. Replacing the inequality with an equality for each of the z'th rows of A, the 
plane of the zth face of the polyhedron is then expressed by 

a n x + a i2 y + a i3 z -b { =0 (12) 

In this way, the region bounded by the intersection of all F infinite planes defined by Equation (11) 
defines the interior region and surface of the polyhedron. Some specific examples are shown below in 
Equations (13) to (16) and Figures 8 to 11. 



Figure 8. — Four-faced polyhedron (tetrahedron) 
particle shape. 
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Figure 10. — Six-faced polyhedron (cube) particle shape. 
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Simulated particle geometries based on polyhedra are a radical departure from the previous shape 
methodologies. General polyhedra can easily be defined with numerous faces, edges, and vertices, as well 
as asymmetries to closely agree with actual granular particle shapes. Traditionally, polyhedra shapes have 
been difficult to implement in DEM software, since the tradeoff between particle shape complexity, 
number of particles, memory, and CPU speed has driven most commercial DEM software towards 
utilizing nonpolyhedra particle shapes. Advances in computer hardware in recent years, however, has 
made utilization of polyhedra particles a practical strategy. The DEM software BLOKS3D developed at 
the University of Illinois by Zhao (Ref. 14), is an example of this strategy. Simulations using 800 
particles composed of tetrahedrons (four vertices), pyramids (five vertices), cubes (eight vertices), and 
rhombic-dodecahedrons (14 vertices) have been reported by Zhao. 

2.8 Contact Detection and Force Model 

The most difficult parts of DEM from a CPU resources perspective are contact detection and contact 
force model implementation (Ref. 15). The complexity of a contact detection algorithm increases rapidly 
with particle shape complexity. Particle shapes that are simple spheres require the least complex detection 
strategy. Poly-ellipsoids, super-ellipsoids (also called superquadrics), asymmetric superquadrics, and 
general polyhedra require more sophisticated contact detection algorithms. 

On the two extremes of contact detection algorithms, there are (1) continuous math algorithms based 
on analytic models of particle shape and (2) discrete math algorithms, which use gridded representations 
of particles. Numerous combinations of analytical and gridded algorithms, along with ad hoc algorithms, 
led to many possible approaches to contact detection. The details of the approach used may be intimately 
dependent on other unrelated aspects of the DEM code. 
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One method of determining the contact force direction in a continuum domain for arbitrarily shaped 
particles, which is loosely based on the method of Hogue (Ref. 16), makes use of the equations of the 
continuous particle surfaces. In this approach, two particle surfaces, described by f (x,y) = 0 and 
g(x, y) = 0 , are said to overlap when a set of coordinates can be found that simultaneously satisfy both 
equalities. For example, the two-dimensional case of Equation (3) defines an ellipse, so that 


f(u, v) = 1 - 


M 

2 

f v l 

\ a \ ) 


V*1 V 


g(«,v) = l- 


r U —Uq ^ 

2 

0 

1 

a 2 J 


l b 2 ) 


(17) 


where u and v are local coordinates of the particle coordinate system. Since the particles are rotated and 
translated relative to the origin of the global coordinate system: 
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(18b) 


then 


f(x,y) = 1 - 
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^xcos^ + ysin0-w o ^ 

2 
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V t>2 2 


(19a) 


(19b) 


For this example with a\ = 0.5, b\ = 1.1, a 2 = 0.4, b 2 = 1.5, u 0 = 0.3, v 0 = 1.8, and c|) = -30°, the 
solution of the simultaneous equations f(x, y) = 0 and g(x, y) = 0 are 




f 0.395809) 



r 

UJ 


v 0.6721 19 j 


Vy 2 ) 

V 


0.0965897 


(20a) 


This case is depicted in Figure 12 where a line segment connects the two solutions. The contact force 
direction angle is the perpendicular to the slope of the line 

\ 

71 


0 = tan 1 


V*2 -*1 J 


+ ■ 


(20b) 


and is equal to 10.08° in the above case. 

In general, there will be two point solutions in two dimensions and a set of solutions comprising a 
closed loop in three dimensions. The two-dimensional case shown in Figure 12 defines a line connecting 
the two points of intersection. A normal vector to this line, located at the midpoint of the line, is one 
method of defining an effective contact force direction as an input prelude to the contact force model. The 
contact force model uses the overlap area or volume information to determine the magnitude of contact 
repulsion. 
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Figure 12. — Contact detection of two particles 
described by analytical expressions. The contact 
normal direction is 0 = 10.08°. 


In the three-dimensional case, the force direction of contacting particles is easily defined for spherical 
particles. Most DEM software, such as PFC3D (Ref. 5b) by Itasca, defines the force direction, in the case 
of simple spherical particles, as along the axis connecting contacting particle centers. In the case of more 
complex particles composed of clusters of overlapping spheres of different sizes, the force direction can 
still be modeled as along the axis connecting sphere centers of the two local spheres involved in the 
contact, assuming that the deformation effect due to contact can be isolated to only the two contacting 
spheres. In this case, the torque about the center of mass of the cluster must also be included in the 
equations of motion. 

The details of force direction are intimately coupled to the details of the particle shape model. For 
example, in the shape model used by Hogue (Ref. 16), particles are made up of triangular surfaces and 
nodes, where the points of the triangles are nodes. When two particles contact, a node of one will usually 
penetrate the triangular surface of another. The force direction is defined as the normal to the triangular 
surface. Other possibilities include node-to-node and triangle-to-triangle contact, which are treated as 
special cases. 

A common characteristic of most approaches to modeling contact force direction is the explicit or 
implicit need to know the normal to the surface of the particle at the point of contact. The gradient of the 
particle surface is a convenient means to quantify the surface normal of the continuous particle shape. For 
example, going back to the ellipsoid of Equation (3), the normal direction at any point on the surface is 
found from the gradient: 
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V f(x, y, z) = 


( 21 ) 
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The difficulty of applying Equation (21) directly to the general case of arbitrary shaped particles is 
that it must be decided at what point on the contacting surface and associated normal direction is most 
relevant to the final contact force direction. A brute force method could involve integrating Equation (21) 
over the overlapping surface areas of both contacting particles to get an average normal direction, and 
thus a final force direction. 

On the other extreme of contact detection methods, the full discrete approach to contact detection can 
be performed using any particle shape since once the shape has been defined, the particle location is 
simply designated by a cluster of pixels. As shown in Figure 13, contact detection is the simultaneous 
enabling of pixel bits in the quantized space of the overlapping particles. The advantage of full discrete 
contact detection is that there is no penalty due to particle shape. However, the obvious disadvantage is 
resolution of the particle states (position and velocity) in that discrete space. 

The discrete version of the contact algorithm involves rectangular gridding of the particles so that any 
overlap can be detected and described by simultaneous occupancy of adjacent particles in one or more 
grids (or pixels for the two-dimensional case). Once a contact condition is detected by the overlap of one 
or more pixels associated with a particle pair, the direction of the contact normal can be computed by 
flagging the set of overlapping pixels. As is true for most mathematical models that substitute a simple 
solution for a more complex mathematical description, numerous approaches are possible. Below we 
derive (as an example) a method for calculating contact force direction that can be implemented as a 
computer algorithm. However we do not imply that this method is necessarily a method of choice by any 
specific DEM software vendor. 

The set of M-occupied pixels X* can be transformed to a principal component direction by using the 
Hotelling transform, and then by taking the perpendicular direction, the contact normal direction can be 
defined. The procedure is described below. 

First, the center of mass m x of the set of overlapping pixels is calculated: 

, M 

= <22) 


Next, the centered covariance matrix is computed: 


, M 

C = — EX -Xf -m - mf 


i=i 


(23) 


The eigenvalues of C x define the principal axis directions. Adding 90° to the angle that defines the 
first principal axis defines the normal direction, which is the contact force direction. 

Several commercial off-the-shelf (COTS) DEM software packages are available. In addition, several 
custom software packages have been developed, some of which may become COTS products in the 
future. In the final analysis, the choice of particle shape is predetermined by the details of specific DEM 
software and the particle shapes it supports. Some of these are listed in Table I. 
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Figure 13(a). — Contact detection of two particles Figure 13(b). — Contact detection of two particles 

described by discrete particle representations described by discrete particle representations 

with fine spatial grid. with coarse spatial grid. 


TABLE I.— PARTICLE SHAPE CAPABILITIES OF SOME DEM SOFTWARE 


Software 

Supported particle shapes 

PFC2D/3D (Ref. 5) 

Spheres, rigid sphere clusters: overlapping and nonoverlapping 

EDEM (Ref. 6) 

Spheres, sphere clusters: overlapping and nonoverlapping 

BLOKS3D (Ref. 16) 

Polyhedra 

ROCKS3D 

Polyhedra 

MIMES (Ref. 17) 

Superquadric shapes 

YADE (Ref. 18) 

Spheres 

Polyhedra (work in progress) 


3.0 Analytical Representations of Particle Size Distribution 

The data-fitting procedure outlined in Bagnold (Ref. 20) provides valuable insight into the nature of 
the desert sand particle size distribution (PSD). The technique used by Bagnold, which was appropriate 
for the time period, assumes that the tools of choice are a pencil, graph paper, and a ruler. Fortunately, 
current computer software tools are beginning to provide a reliable means to fit arbitrary nonlinear 
functions to a data set. Taking this approach, it makes sense to define the PSD as a single formula, 
covering the entire particle size range. 1 


^hen dealing with particle size histograms and the numerous forms of size distributions and their moments, 
carefully keeping track of units reduces confusion and the risk of error. Also, it should be noted that histograms are 
summed, while distributions are integrated. 
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3.1 Size Distribution by Mass 

The fractional mass of particles in a differential size increment d D, based on Bagnold’s work is 


N(D)dD = - ± d D (24a) 

((A>/f>r /T +(A)/£>r /Y ) y 

where the small slope s is always positive and coarse slope c is always negative. The PSD of 
Equation (24a) is normalized by setting the integral of N(D ), for sizes from 0 to °o ? equal to 1 : 


J N(D)dD = 1 (24b) 

0 

Note that in the case of Bagnold’s data, the integral is set equal to 100 percent. 

For the small-diameter case, D « Z) 0 , so that Equation (24a) reduces to 


N(D) = N 0 (D/D 0 ) S 


(25) 


If Equation (25) is plotted on a log-log plot, as described by Bagnold, a straight line graph of slope 5 
results in 


y = log TV = log7Vo +sk>g J D- l sTogD 0 
= y o + sR - sR 0 


(26) 


where R = logD . 

In the large diameter case, D » D 0 , so that Equation (24) reduces to 


N(D) = N 0 (D/D 0 ) C 


(27) 


If Equation (27) is plotted on a log-log plot, a straight-line graph of slope c results, which is 
equivalent to 


y = logTV = log7V 0 +clogD-clogD 0 
= y 0 +cR-cR 0 

Recall that c has been defined as negative value according to Bagnold. 

The exponent parameter y in Equation (24) is introduced as a means to control the transition rate 
between the small-diameter positive slope regime and the large-diameter negative slope regime. 
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Grain diameter, D [mm] 


Figure 14. — Particle size distribution fit to desert sand using 
Equation (24a). 


3.3 Particle Size Distribution 

Equation (24) describes a size distribution based on particle mass. This formula is an empirical fit to 
observed grain size distributions of desert sand (Fig. 14) based on measurements using a multilevel sieve. 
However, the fundamental size distribution of primary interest for an assemblage of particles is P(D ), 
which may be expressed in units of [m -3 mm -1 ], where P(D ) d D is the number of particles per unit volume 
of space within a differential size increment dD. The total number of particles of all sizes in a sample 
volume Ks-is then 


N v 



(29) 


Assuming the packing density is homogeneous and repeatable, the total mass of the assemblage of 
particles in a volume V s is 


M v = \npV s | D 3 P(D)dD (30) 

0 

where p is the particle grain density, assuming all grains have equal density and are spherical. 

Equation (30) can be rewritten in terms of the volume density of the assemblage of particles by dividing 
both sides by V s : 


p v =}npfD 3 P(D)dD 

0 


(31) 
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Now the mass size distribution described by Equation (24) can be expressed in terms of P(D ): 


N(D) = 


jnpD^PjD) 

}7tp Jz) 3 P(D)dZ) 

0 


= ±kD 3 P{D)-^ 
Pv 


( 32 ) 


Equation (32) can now be inverted, thus providing a means to estimate the particle size distribution 
P(D ) from the measured mass distribution N(D ): 

P(D) = t ^-—N(D) 

I nD 3 p 

(33 

_ P v No 

6 nI)3 P ((Z) 0 / D)/y +(D/ D 0 )~/v f 

where the parameters No, 5 , c , y, p, and p^are determined by measurement. Note that the particle packing 
ratio is less than one (p v /p < 1) , because of gaps between particles of an assemblage. 


3.4 FPA Data of JSC-la 

Instruments that measure and count individual particles typically output a histogram count of particles 
in bins of equal intervals for the property of interest. This is the case with the fine particle analyzer (FPA) 
used to measure a sample of JSC-la at the NASA Kennedy Space Center Granular Mechanics and 
Surface Systems Laboratory. An example of particle size counts from the first nine size bin intervals of an 
FPA measurement are shown in Table II. Each bin interval is defined by a AD = 1/ 256 [mm], shown in 
the first column. The second column shows counts C k of particles that fall into the size bin centered at the 
equivalent spherical diameter D k where k is the bin index. 

To convert from a histogram of counts to a size distribution n{D) in units of [mm 1 ], the counts are 
divided by the bin interval, AD : 


n (D k ) ~ n k = %- [mm- 1 ] (34) 

AD 

The third column of Table II shows n k , , the corresponding values of the number of particles in bin k , 
determined from the C k and AD = 1/ 256 [mm], using Equation (34). The n k are plotted in Figure 15, 
along with an empirical fit based on a sum of exponentials: 

n(D) = + N 2 e~ D//jyi [mm 1 ] (35) 
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TABLE II.— FPA a PARTICLE COUNT DATA FOR JSC- la 


Particle size, 
D k 
[mm] 

Particle count, 

Q 

Count per size interval, 
n k = C/JAD 
[mm 1 ] 

0.003906 

1232155 

315431680 

.007813 

629101 

161049856 

.011719 

676994 

173310464 

0.015625 

659910 

168936960 

.019531 

602111 

154140416 

.023438 

546315 

139856640 

0.027344 

496099 

127001344 

.03125 

448268 

114756608 

.035156 

403354 

103258624 


a FPA refers to fine particle analyzer at NASA Kennedy 
Space Center. 


TABLE III.— FPA a MASS DATA FOR JSC- la 


Particle size 
Dk , 

[mm] 

TL ratio, 
thickness/length 

Mass per interval, 
[gmirf 1 ] 

0.003906 

0.30649069 

0.005223 

.007813 

.51233492 

.041281 

.011719 

.54465904 

.162795 

0.015625 

0.56233665 

0.392799 

.019531 

.57190647 

.716264 

.023438 

.57940036 

1.143148 

0.027344 

0.58504106 

1.670413 

.03125 

.58960083 

2.277133 

.035156 

.59410737 

2.948024 


a FPA refers to fine particle analyzer at NASA Kennedy 


Space Center. 



Figure 15. — Fine particle analyzer (FPA) size distribution 
of JSC-1 a, fitted to sum of exponentials. 


Particle mass data corresponding to the counts and diameters of Table II are shown in Table III. In the 
second column, L is defined as the longest chord across the particle, where a chord is defined as any line 
segment that crosses from one part of a particle's perimeter to another part of the same perimeter; and T is 
defined as the longest chord of all chords that are perpendicular to L. Typically in this context L denotes 
particle length and T, thickness. 

The size distribution is converted to mass distribution by multiplying by particle volume and particle 
density: 


m(D) = ±nD 3 p Rtl(1 + Rt, ) n(D) [g mm -1 ] (36) 

where the particle density p = 2.65 [g em -3 ]. The geometric factor R TL (1 + R TL ) / 2 in Equation (36) 
improves the particle volume estimate using a geometric correction factor to account for the nonspherical 
nature of JSC-la. The TL ratio, R TL shown in Table III, is an output of the FPA data processing. Note 
that as R tl approaches 1, the geometric correction factor also approaches 1. 


NAS A/TM— 20 1 0-2 1 6257 


18 


Using the empirical n(D) function from Equation (35) and the equation for mass shown above, the 
FPA particle mass data is plotted in Figure 16. 

Going back to Bagnold’s method, Equation (24), the particle mass data from Table III can be fitted as 
shown in Figure 17. In this case the parameter N 0 is expressed so that P(D ) from Equation (29) is in units 
of [g-mnT 1 ]: 


/ 77 (D) = 


Nc\ 


{DqId)7i+{d 0 



[gmm ‘] 


(37) 


A small-slope parameter s and coarse-slope parameter c is determine by this fit. The FPA particle 
count size distribution n k is plotted in Figure 18 along with the Bagnold’s fit, found by inverting 
Equation (36) with m{D) from Equation (37): 


n(D)= Wp 


1 _ n 3^ R TL (1 + R TL ) 


m(D) 

No 


[ mm " 1 ] 


1 D 3 n r tl(} + r tl) 

g KL> p 2 


({Do!D) sh + {DIDo)~ cl ^ 


For comparison, a lognormal fit to the mass size distribution is given by 

M ) 2 

m$ e v uy 


V) 


/ 77 (D) = 


Dq ^a(D/D 0 ) 


[g mm *] 


(38) 


(39) 



D [mm] 

Figure 16. — Fine particle analyzer (FPA) mass size distribution of JSC-1 a, fitted to sum of 
exponentials, Equation (35). 
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m(D) [g • mm 



D [mm] 

Figure 17. — Fine particle analyzer (FPA) mass size distribution of JSC-1 a, fitted to Bagnold’s model. 



Figure 18. — Fine particle analyzer (FPA) size distribution of 
JSC-1 a, fitted to Bagnold’s model using Equation (38). 
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D [mm] 

Figure 19. — Fine particle analyzer (FPA) mass size distribution of JSC-1 a, fitted to lognormal distribution. 


The fit of Equation (39) to JSC- la is shown in Figure 19. 


4.0 Concluding Remarks 

This report discusses topics that are important to the modeling of the interaction and dynamics of 
granular particles, and in particular, lunar stimulants. In Section 2.0, discrete element method (DEM) 
particle shapes were presented, where the choice of particle shape is motivated by many opposing 
tradeoffs, such as 

(1) The number of parameters describing the simulated particle surface (thus defining the interior 
versus exterior region). The goal is to attempt to closely match the shape of a real particle. 

(2) The total number of particles that can be incorporated in a simulation. 

(3) Computer processing resources needed: memory and total simulation time. 

(4) The difficulty versus ease of the contact detection algorithm and subsequent implementation of a 
force model. 

Most commercial and custom DEM software supports some type of complex particle shape in order 
to introduce asymmetry into the particle shape, so as to better simulate real granular particles interactions 
and behavior. 

In Section 3.0, a complete particle size distribution (PSD) formula is shown to fit the example desert 
sand data from Bagnold, using a coarse-grade exponent c and small-grade exponent 5. An additional 
empirical parameter y has been introduced to connect the two log-log linear fits in the small- and coarse- 
grade particle size regions. The key finding in this section is that Bagnold's model is a fair representation 


NAS A/TM— 20 1 0-2 1 6257 


21 


of the PSD except for particles below 20 pm, where it diverges somewhat. Also, the double exponential 
formulation is a fair representation of the size distribution, but the lognormal formulation is not very 
successful at matching the JSC- la lunar stimulant used in this evaluation. Particle size data of JSC- la 
obtained from a fine particle analyzer (FPA) at NASA Kennedy Space Center are also fitted to the PSD in 
order to extract equivalent small- and coarse-grade exponents. Note that the FPA data are uncorrelated 
quantitatively with other accepted PSD measuring tools, so there may be a bias in data that affects the 
quality of the PSD fit. 
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Appendix A. — Microscope Image, Shadowgraphs, and 
SEM Spectrum of JSC-la 

JSC- la is a lunar regolith simulant that was developed by the NASA Johnson Space Center. It is a 
basaltic ash with a high glass content. Figure 20 is an optical microscope image of the larger particles in a 
sample of JSC-la. Figure 21 presents shadowgraphs of JSC-la stimulant showing the variety of shapes 
and sizes in the mixture. Figure 22 is an elemental analysis spectrum of a sample of JSC-la using a Joel 
75000F scanning electron microscope (SEM) with energy-dispersive x-ray spectroscopy (EDS). Figure 23 
contains elemental maps of the sample in Figure 22, showing K-a emission lines for the detected 
elements (K lines result when an electron transitions to the innermost K shell from a 2p orbital of the L 
shell). 



Figure 20. — Microscope image of >600-pm particles of JCS-la. 



Figure 21 . — Shadowgraphs of JCS-la. 
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Full scale counts: 4857 Extracted Spectrum Cursor: 1.619 keV 


261 Counts 



klm - 56 - Bci keV 

Fri Jul 11 14:09:19 200S 
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Al K 
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0.9S 
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0.02 

Ti K 
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0.17 
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0.01 

FeK 

3.S5 
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1.05 

+/- 

0.03 

Cu K 

1.93 
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0.46 
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0.05 
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0.03 
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100.00 
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Figure 22. — Joel 75000F scanning electron microscopy (SEM) spectrum of JSC-1 a. 
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Figure 23. — Energy-dispersive x-ray spectroscopy maps, showing spatial mapping of elements corresponding to 
scanning electron microscopy (SEM) spectrum of Figure 22. 
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Appendix B. — Summary of Particle Shape Considerations 

Discrete element method (DEM) simulation software considers a granular material as being made up 
of a size distribution of idealized particle shapes. These shapes may be simple two-dimensional 
nonoverlapping or overlapping discs or ellipses, three-dimensional spheres or ellipsoids that overlap or do 
not overlap, or complex shapes such as polyhedrons or superquadrics. Table IV compares some of the 
pros and cons of these idealized particle types. 


TABLE IV.— COMPARISON OF DISCRETE ELEMENT METHOD (DEM) PARTICLE SHAPE TYPES 


Particle shape 

Pros 

Cons 

Comments and possible 
solutions to cons 

Incorporated in DEM 
software 

Sphere 

Simplest to 
implement; requires 
the minimum CPU a 
resources, allowing a 
maximum number of 
simulated particles. 

Particle behavior is 
unrealistic because 
of simplified shape. 

Implements more complex 
force and friction models 
to match real particle 
behavior. 

PFC2D/PFC3D (Ref. 5) 
EDEM (Ref. 6) 

Yade (Ref. 19) 

Rigid 

nonoverlapping 
sphere cluster 

Relatively simple to 
implement. 
Approximates real 
particle shapes. 

In order to approximate 
a real particle shape, 
many spheres of 
various sizes are 
required, significantly 
increasing CPU 
resources. 

Decreases number of free- 
state parameters to only 
that needed to simulate a 
single rigid body. 

PFC2D/PFC3D (Ref. 5) 
EDEM (Ref. 6) 

Rigid 

overlapping 
sphere cluster 

Adds ability to better 
match real particle 
shapes. 

Again, in order to 
approximate a real 
particle shape, many 
spheres of various sizes 
are required. Some 
DEM codes (e.g., 
EDEM) can correctly 
account for total cluster 
mass density. 

Specifies a variable sphere 
mass density so that total 
cluster is the correct value. 
This increases complexity 
of sphere density 
specification. 

PFC2D (Ref. 5a) 

Ellipsoid 

Better matches real 
particle shapes. 

Contact force 
detection has some 
similarities to that of a 
sphere. CPU resources 
required are 
comparable to that of a 
sphere. 

Contact force detection 
is more complicated 
than for a sphere. The 
ability to match real 
particles is much better 
than a sphere, but is 
still limited. Contact 
detection is more 
difficult than for 
sphere. 

Ellipsoid is comparable to 
rigid sphere cluster as far 
as tradeoffs between CPU 
resources, number of 
possible particles, and 
similarity of simulated to 
real particle shape. 

Wang (Ref. 21) 

Superquadric 

More simulated 
particle parameters are 
provided in a closed 
form description, 
which result in more 
shapes. All of the 
particle parameters 
can be computed 
analytically. 

Similar difficulties 
encountered as with 
ellipsoids. When 
concave exponent 
values are used ( n < 2), 
additional difficulties in 
contact detection and 
force model are 
present. 

Because purely analytic 
methods used for contact 
detection and force model 
can get messy, a gridded 
spatial model of particle 
and particle interactions 
may be a better choice. 

MIMES (Ref. 17) 

Asymmetric 

superquadric 

Similar to 
superquadric, but 
many more parameters 
can be added, 
providing a better 
match to real particle 
shapes. 

All problems described 
for superquadric 
particle shape are 
compounded by the 
inclusion of additional 
parameters. 

Even though an analytical 
approach to contact 
detection and force model 
is possible, it may not be 
practical in this case. A 
gridded numerical strategy 
may be the best approach. 

None found 
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Particle shape 

Pros 

Cons 

Comments and possible 
solutions to cons 

Incorporated in DEM 
software 

Eight-quadrant 

superquadric 

A subset of 
asymmetric 
superquadric particle 
type, which provides a 
method to easily 
guarantee no surface 
discontinuities. 

Problems are similar to 
the asymmetric 
superquadric. 

Gridded representation of 
particle and particle 
interactions is probably the 
most practical approach. 

None found 

Poly-ellipsoid 

A subset of eight- 
quadrant superquadric 
particle type, which 
guarantees no surface 
gradient 

discontinuities (i.e., no 
sharp ridges or edges). 

Problems with contact 
detection are more 
severe than with the 
ellipsoid, but less 
severe than with the 
general superquadric. 

A gridded approach may 
be best, but other pseudo 
analytical approaches have 
been developed which 
claim good success. 

Peters (Ref. 13) 

Polyhedra 

The polyhedron-based 
particle shape is the 
most general shape of 
those discussed in this 
summary. 

Analytical closed form 
contact detection 
methods are very 
difficult except in 
specific simple cases. 
Contact force 
determination on edges 
is also difficult and 
ambiguous. 

A gridded approach is 
almost a necessity. Sharp 
edges can be eliminated by 
performing a spatial 
convolution of the particle 
surface or volume with a 
sphere (Minkowski sum, 
Ref. 22). 

Nezami (Ref. 15) 


a CPU refers to central processing unit. 
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Appendix C. — Eight-Quadrant Superquadric Particle Shape Examples 

Using Equation (8), the 12 particle shape parameters are randomly modified to demonstrate the 
resulting diversity of particle shapes. The figures in this appendix were created using Mathematica 
version 7.0. 


x(0, (])) = < 
7(e,4>) = < 

z(0, ())) = < 


^(sinGcoscj )) 2 ^ 1 
a 2 |sin0cos())| 2/i:2 
b\ (sin©sin(})) 2/l1 
Z> 2 |sin 0sin c))| 2/Z ' 2 
C| (cos0) 2/ W| 
c 2 |cos0| 2/A/2 


for 

for -|<<| )<4p 
for 0 < (|) < 7i 
for n < (|) < 271 
for O<0<-| 
for j < 0 < 71 


O<0<7t 


O<0<71 
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